
import indexing.featureindex
import sys
import coverage.readsperbasefromsummary
from genome.genomereader import GenomeReader

readers = {}

def getReadsAsList(reads, genomic):
    if "/" in reads:
        return reads.split("/")
    if genomic=="X":
        ix=0
    else:
        ix = list("ACGT").index(genomic)
    ret = ["0", "0", "0", "0"]
    ret[ix] = reads
    return ret

def getCoverage(index, chr, pos):
    features = list(index.getFeatures((chr, pos, pos)))
    if not features:
        # ['T\t34\t0.853571\t27\tT\t0/0/0/28\t0/0/0/24\n']
        content = "X\t0\t0.00\t0\tX\t0\t0\n"
    else:
        content = features[0].content
    (genomic, numReads, weightedReads, numFragments, isSnp, readsByBase, fragmentsByBase) = tuple(content.strip().split("\t"))
    readsByBase = "\t".join(getReadsAsList(readsByBase, genomic))
    fragmentsByBase = "\t".join(getReadsAsList(fragmentsByBase, genomic))
    return "%s\t%s\t%s\t%s" % (numReads, numFragments, readsByBase, fragmentsByBase)

def getComplexCoverage(indexFilename, chr, pos, newBase, filtered=True):
    if not readers.has_key(chr):
        readers[chr] = GenomeReader(chr)
    reader = readers[chr]
    features = list(coverage.readsperbasefromsummary.process((chr, '?', pos, pos), indexFilename))
    if not features:
        feature = coverage.readsperbasefromsummary.ReadsPerBase()
    else:
        (p, feature) = features[0]
    numReads = feature.getReads(filtered)
    numFragments = feature.getFragments(filtered)
    reads = "\t".join(map(str, feature.getReadCounts(filtered)))
    fragments = "\t".join(map(str, feature.getFragmentCounts(filtered)))
    homopolymerLeft = 0
    homopolymerRight = 0
    seq = newBase
    while homopolymerLeft < 100 and reader[pos-homopolymerLeft-1] == newBase:
        homopolymerLeft += 1
        seq = newBase.lower() + seq
    while homopolymerRight < 100 and reader[pos+homopolymerRight+1] == newBase:
        homopolymerRight += 1
        seq = seq + newBase.lower()
    for i in xrange(5):
        seq = reader[pos-homopolymerLeft-1-i].lower() + seq + reader[pos+homopolymerRight+1+i].lower()
    return "%ld\t%ld\t%s\t%s\t%dL\t%dR\t%s\t%d" % (numReads, numFragments, reads, fragments, homopolymerLeft, homopolymerRight, seq, 1+homopolymerLeft+homopolymerRight)

if __name__=='__main__':
    useCoverageSummary = False
    normalsamp = sys.argv[1]
    if useCoverageSummary:
        ixname = "/archive10/bowtie/Full_Length_Results/Sample_Parameters/%s/Current/%s.coverage_summary.INDEX" % (normalsamp, normalsamp)
        print >> sys.stderr, "# Using coverage summary file", ixname
    else:
        if "/" in normalsamp:
            ixname = normalsamp
        else:
            ixname = "/archive10/bowtie/Full_Length_Results/Sample_Parameters/%s/Current/%s.bowtie_sorted_summary.INDEX" % (normalsamp, normalsamp)
        print >> sys.stderr, "# Using sorted summary file", ixname
    ix = indexing.featureindex.FeatureIndex(ixname)
    # print >> sys.stderr, "# WARNING: filtering for dbsnp-yes, frac >= 0.25, var frags >= 8"
    # print >> sys.stderr, "# WARNING: filtering for dbsnp-no, frac>=0.1"
    # print >> sys.stderr, "# WARNING: filtering for frac>=0.1; no filter on dbsnp yes/no"
    print >> sys.stderr, "# WARNING: no filter on frac or dbsnp yes/no"
    for line in open(sys.argv[2]):
        cols = line.strip().split("\t")
        # if cols[16]!="dbsnp-no" or float(cols[4])<0.1 or (cols[18]!="non-coding-3-prime-utr" and cols[18]!="non-syn"):
        # if cols[16]!="dbsnp-no" or float(cols[4])<0.1:
        #    continue
        # if cols[16]!="dbsnp-yes" or float(cols[4])<0.25 or int(cols[3])<8:
        #    continue
        #if float(cols[4])<0.1:
        #    continue
        snpId = cols[1]
        fields = snpId.split("_")
        (chr, pos, oldBase, newBase) = (fields[0], long(fields[1]), fields[2], fields[3])
        if useCoverageSummary:
            ncov = getCoverage(ix, chr, pos)  # use this to read from coverage map file
            print "%s\t%s" % (ncov, line.strip())
        else:
            ncov = getComplexCoverage(ix, chr, pos, newBase) # use this to read directly from bowtie sorted summary file
            cols = ncov.split("\t")
            print "%s\t%s\t%s" % ("\t".join(cols[:-4]), line.strip(), "\t".join(cols[-4:]))

